Spin-Dynamical Analysis of Supercell Spin Configurations 



NO 
O 
O 
(N 

in 



5h 

43 



Andrew L. Goodwin, 1 Martin T. Dove, 1 Matthew G. Tucker, 2 and David A. Keen 2,3 

1 Department of Earth Sciences, Cambridge University, Downing Street, Cambridge CB2 3EQ, U.K. 
ISIS Facility, Rutherford Appleton Laboratory, Chilton, Didcot, Oxfordshire 0X11 OQX, U.K. 
1 Department of Physics, Oxford University, Clarendon Laboratory, Parks Road, Oxford 0X1 3PU, U.K. 

(Dated: February 6, 2008) 

A model-independent approach capable of extracting spin-wave frequencies and displacement vec- 
tors from ensembles of supercell spin configurations is presented. The method is appropriate for 
those systems whose spin-dynamical motion is well characterised by small-amplitude fluctuations 
that give harmonic spin waves. The generalised spin coordinate matrix — a quantity that may be cal- 
culated from the observed spin orientations in an ensemble of spin configurations — is introduced and 
its relationship to the spin-dynamical matrix established. Its eigenvalues are subsequently shown 
to be related to the spin-wave mode frequencies, allowing the extraction of spin-wave dispersion 
curves from configurational ensembles. Finally, a quantum-mechanical derivation of the same re- 
sults is given, and the method applied as a case study to spin Monte-Carlo configurations of a 3D 
Heisenberg ferromagnet. 
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I. INTRODUCTION 

"Atomistic" simulations of spin interactions in mag- 
netic systems provide a useful mechanism of understand- 
ing both structural and dynamical properties of this im- 
portant class of materials. The reverse Monte-Carlo 
(RMC) refinement method is a relevant example: a data 
driven approach, RMC attempts to account for observed 
experimental data in terms of ergodic assemblies of atom- 
istic and spin configurations iSi 3 . There are a number 
of established methods of analysing such configurations 
so as to calculate important thermodynamic quantities; 4 
however, to the best of our knowledge a method of cal- 
culating spin- wave dispersion curves in this way has not 
yet been reported. In this paper, we describe one such 
method. An important feature of this approach is that 
spin-dynamical quantities might be extracted from a vari- 
ety of experimental methods, if these are used to generate 
the configurations (e.g. by using RMC). Indeed, we have 
recently applied this approach to reverse Monte-Carlo 
configurations, allowing extraction of spin-wave disper- 
sion curves from neutron diffraction data£ Our method 
is appropriate for those systems whose spin-dynamical 
motion is characterised by small-amplitude fluctuations 
that give harmonic spin waves. We require harmonic- 
ity only at constant temperature; some degree of anhar- 
monic behaviour is allowed in the sense that it is possible 
to observe changes in spin-wave energies across analyses 
corresponding to different sample temperatures. 

Our paper is arranged as follows. In section II, we 
first briefly review the well-established concept of the 
spin dynamical matrix, re-casting some of the key equa- 
tions in an appropriate form for our supercell configu- 
rational analysis. We then introduce the spin coordi- 
nate matrix, which is assembled from the normal mode 
spin coordinates. We proceed by using this matrix to es- 
tablish a relationship between the observed orientations 
in real-space spin configurations and the spin-wave fre- 



quencies. Furthermore, we show explicitly how one can 
then extract spin-wave dispersion curves and spin-wave 
mode assignments from ensembles of spin configurations, 
without any prior knowledge of the exchange constants. 
Finally, we present a derivation of the key results from 
a quantum-mechanical perspective. In section IIIII we 
illustrate the capabilities of the approach through the 
analysis of configurations generated using a simple spin 
Monte-Carlo simulation. We show that the method yields 
the expected spin-wave dispersion relation. Finally, sec- 
tion IV discusses the implementation and some possible 
extensions of the method. 



II. THEORY 



Semi-classical derivation of the "spin coordinate 
matrix" 



Our theoretical approach builds on the well-established 
semi-classical the ory o f spin dynamics (as described, for 
example in Refs.lf H7l8l) . Throughout our analysis we will 
consider only systems with a single spin-alignment axis; 
however, an extension to multiple-axis systems would be 
straightforward. Our starting point is the spin-dynamical 
matrix A(k), which stores dynamical information in 
terms of the Heisenberg exchange integrals J. Its rows 
and columns are indexed by the spin types j, and the 
corresponding elements are given explicitly by the ex- 
pression:' 
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Here, the Sj and r{j£) are the spin quantum numbers and 
average positions, respectively, of each spin j in each unit 
cell I. The eigenvalues of the spin-dynamical matrix are 
the spin- wave frequencies w(k, v), and the right eigenvec- 
tors describe the normal mode spin-wave displacements. 6 

We now construct the matrix S(k), which we will term 
the "spin coordinate matrix" . It is built from the indi- 
vidual spin orientations produced by the set of spin- wave 
modes acting on a spin configuration. Two important 
results will be derived: first, that the spin coordinate 
matrix is related to the spin-dynamical matrix by the 
appealing relationship A(k) • E(k) = k B T (in the clas- 
sical limit); second, that the eigenvalues of S(k) are the 
quantities k^T/cufe, v). The latter property is the key 
result of this paper, as the spin coordinate matrix can be 
calculated directly from spin configurations, and conse- 
quently the spin-wave frequencies can be determined via 
matrix diagonalisation. 

We begin by introducing the magnon variables 
r(j, k, t), r*(j, k, t). These are defined in terms of the 
standard spin oscillators a + and 

r(j,M) = W^^a+(^,t)exp[ik.r(j£)], 

e 

r*(?,M) = W^-^CT-(^,t)e X p[-ikT(^)]. 



These conjugate variables may be calculated directly for 
each spin configuration (indexed by the variable t). In 
particular, they do not rely on any a priori knowledge 
of the exchange integrals J, nor the number and type of 
significant interactions for each spin. We proceed by as- 
sembling the magnon variables into the column vectors 
r(k, t), T*(k, t), from which the t-averaged spin coordi- 
nate matrix is formed: 

£(k) = (x*(k).T T (k)). 

The elements of S(k) are given by 



the product matrix are given explicitly by: 

[A(k) • S(k)] iJ = £ S 3 S r jf JJ '\ [(*-(jO)v + W) 
j>,e< ^ ' 

(3) 

The off-diagonal entries [A(k) • E(k)]jj (i ^ j) con- 
tain terms whose exchange integrals J{i,j') and spin dis- 
placement correlations (<J~ (j')a + (j)) do not share like 
terms; consequently, their time average vanishes. As 
such, A(k) • E(k) is a diagonal matrix whose diagonal 
entries are given by Eq. • 

Within the high-temperature limit, we proceed to show 
that the diagonal entries are equal to the thermal energy 
k^T. This is seen by expanding the spin Hamiltonian 
dot product in terms of the spin variables: 



S 7 ; ■ Si 



S l S 3 [1 



2 W 



Here we have neglected fourth-order terms in as — 
in the limit of small fluctuations — one has \a \ <C 1. 
Substituting this expansion back into the Hamiltonian, 



(4) 

The first term on the right-hand side of this expression 
corresponds to the spin interaction energy of the system 
in the absence of any spin excitations, while the second 
term contains the correction that arises from the popu- 
lation of spin-wave modes. The summation in this spin- 
wave term runs over all non-trivial pairs and so can 
be separated into contributions Hi from each spin i: 



where 



(a+a: 



W (2) 
xexp{ik-[r(/0-r(j£)]}. 

The angular brackets used here represent an average 
taken either over time {e.g., for spin dynamical (SD) sim- 
ulations) or over configurations {e.g., for ergodic ensem- 
bles). In practice, some degree of error will necessarily 
be introduced in this process, by considering finite time 
intervals and/or finitely large configurational ensembles. 
This error diminishes with the square root of the number 
of configurations or timesteps included in the averaging 
process. 

We proceed to investigate the properties of the ma- 
trix S(k) further, first forming its product with the spin- 
dynamical matrix, A(k) • S(k). The diagonal entries of 



Here, each Hamiltonian has been replaced by its time 
average, and so the Hi give the time-averaged thermal 
energy of each spin i; in the high-temperature limit this 
is simply k^T, by the classical equipartition result. A 
more formal quantum-mechanical treatment is given be- 
low; however, the classical result suffices here. By com- 
parison with Eq. J^J i on e has 



[A(k) • E(k)] w = k B T. 



(5) 



This is a key result and gives that the matrix A(k) • E(k) 
is diagonal with diagonal entries k^T . 

One corollary of this result is that the eigenvalues of 
E(k) are the quantities k&T/uj()s., v) for each spin-wave 
mode v. This is an important result because the matrix 
E(k) is constructed from the r(k, t) and r*(k, t), the val- 
ues of which can be determined directly from inspection 
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of the spin displacements in spin configurations. Conse- 
quently, the calculation and diagonalisation of the spin 
coordinate matrix acts as a model-independent method 
of extracting spin- wave frequencies from spin configura- 
tions. 



B. Implementation 

In order to assemble a set of spin-wave dispersion 
curves, one need only calculate S(k) for an appropriate 
range of values of the wave-vector k (corresponding to, 
for example, the particular reciprocal-space directions of 
interest). The maximum "resolution" obtainable relies 
on the number of unit cells represented by each configu- 
ration: a supercell containing n a ,nb,n c unit cells along 
axes a, b, c permits the set of wave- vectors: 

= — a H b -I c z a ,Z6,z c eZ. 

n a n b n c 

Diagonalisation of S(k) at each wave- vector then gives 
the frequencies of the spin-wave modes at that point in 
reciprocal space. But the method yields more in that 
the eigenvectors of S(k) describe the spin-displacement 
patterns associated with each normal mode. In this 
way, the modes may be individually labelled according 
to their particular displacement pattern, and connected 
accordingly from k-point to k-point to form the spin- 
wave dispersion curves. By calculating mode frequencies 
at symmetry-equivalent wave- vectors, it is possible to es- 
timate the error associated with each point on the curve. 
Moreover, the acoustic mode frequency at k = gives 
a first-order estimate of the configurational finite-size ef- 
fects. 

It would be possible to extend this analysis through 
consideration of the implications of lattice symmetry on 
the form of A(k) and S(k). By this we mean that one 
may in principle calculate first the normal mode displace- 
ment vectors using standard symmetry arguments, and 
then use these as a basis in which to express the observed 
spin displacements. In this representation, the spin dis- 
placement matrix would be in block-diagonal form, with 
each block corresponding to a particular irreducible rep- 
resentation of the point group of the magnetic structure. 
Individual blocks could then be separated and treated 
individually, allowing unambiguous identification of the 
symmetry of each of the normal modes. A similar ap- 
proach has been applied elsewhere to the analysis of 
phonon modes from atomistic configurations. 910 



C. Quantum-mechanical derivation of the "spin 
coordinate matrix" 

We now recover the above results from a quantum- 
mechanical analysis, using standard definitions and no- 
tation as described in e.g. Ref.|8|. For convenience, this 
derivation is limited to systems that contain only one 



spin type; however, the formalisms developed are readily 
extended to include multiple spin systems. It is easily 
shown that the spin Hamiltonian can be expressed in 
terms of the standard magnon creation and annihilation 
operators a* , a: 



n 
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This result holds only under quasi-saturation conditions, 
where the spin deviations are smalls Noting the simi- 
larity to Eq. (0} , the first term on the right-hand side of 
this expression represents the energy contribution that 
arises from the spin alignment process itself (i.e., in the 
absence of any spin- wave excitations); the second repre- 
sents the (positive) energy contribution due to thermal 
excitations of the spins. 

We now introduce the Holstein-Primakoff magnon 
variables fii 

b k = ■y=^2& j exp(ik'r J ), 

j 

K = 7^5Z a i ex P(~ ik ' r j)' 

j 

which are related to the creation and annihilation oper- 
ators by reverse Fourier transform: 



-^=^b k exp(-ik-r J ), 
-^=^b k exp(ik-r,-), 



and can be used as a basis in which to express the spin 
Hamiltonian. Some manipulation gives 

H = -J2j(hj)S 2 + J2^KW, (6) 
»W k 

where 

7k = |at^ J ( i >i)[ 1_ex P( _ik-r ii)] 

and Yji = Tj — Yi. Again, it is the second term on the 
right-hand side of Eq. © that arises from the popula- 
tion of spin-wave modes. This expression can in fact 
be described in terms of the spin-dynamical and spin 
coordinate matrices as defined above. The identities 
&j = {S/2) 1 / 2 <j+ and a* = {S/2) 1 / 2 <j- give 

ftb k b k = Ip^ajaj-exptik-r^) 
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Moreover, 



7k 



2S 

hN ^ 



J 



JJ 
it' 



(exp{-ik-[r(iY)-r(^)]}-l) 



= A(k), 

by comparison with Eq. QJ. Consequently, the Hamilto- 
nian can in fact be written in the form 



W=-]Tj(z,j)S 2 + 5>(k).£(k). 



(7) 



The spin-wave energy term in Eq. (JJJ) represents the 
contribution from thermal population of the individual 
modes (treating modes at different wave-vectors sepa- 
rately). The final step in the derivation arises from 
the quantum-mechanical result that the magnon vari- 
ables give the mode occupation number^ the eigenval- 
ues of b£bk are the occupation numbers n(k). Conse- 
quently the eigenvalues of S(k) = fib k bk are the val- 
ues /m(k) ~ &BT/o;(k) (in the high temperature limit). 
This quantum-mechanical analysis also suggests the di- 
agonalised form of S(k) (which we denote as fi(k) and 
A(k) ■ S(k) when the high-temperature approximations 
are not used: 



O(k) = hn(k) 



exp[huj(k)/k B T] - V 



giving in turn 
A(k) • S(k) 



= w(k)n(k) 



huj(k) 



exp[^(k)/fc B T] - 1' 



(8) 



(9) 



We note that in the high-temperature limit, Eq. @ re- 
duces to the classical equipartition result cited earlier. 



III. CASE STUDY: SPIN MONTE-CARLO 
SIMULATION OF A 3D S = § HEISENBERG 
FERROMAGNET 



In order to test our approach we prepared an ensem- 
ble of spin configurations using spin Monte-Carlo sim- 
ulations. Spin Monte-Carlo is an ideal method for a 
test case such as this, as one has complete control over 
the strength and nature of spin interactions, and hence 
the spin dynamical information reflected in its output. 
Consequently the form of the spin dispersion that should 
be recoverable from a configurational ensemble is known 
precisely. For example, the spin- wave dispersion for a 3D 
S = 5 Heisenberg ferromagnet is easily shown to be 

foj(k) = 4SJ[3 — cos(ak x ) — cos(ak y ) — cos(afc z )]. (10) 

Spin Monte-Carlo simulations in which pairs of spins in- 
teract via classical ferromagnetic Heisenberg potentials 



a 
tq 





[oo a 











R 



X 



R 



Wave-vector 



FIG. 1: Spin- wave dispersion curves in a 3D Heisenberg fer- 
romagnet calculated using Eq. 1101 (solid line) , and extracted 
from our spin Monte-Carlo simulations using the method de- 
veloped in the text (open circles — 10 mK; filled circles — IK). 



should then yield configurations that reflect this same 
dispersion behaviour. 

With this in mind, we generated an ensemble of ca 
500 equilibrium Monte-Carlo configurations, each repre- 
senting a 10x10x10 supercell of an idealised primitive 
cubic ferromagnet. The simulation employed a coupling 
constant J — 0.4 meV and was carried out at two temper- 
atures below its paramagnetic phase transition tempera- 
ture T c ~ 3.5 K: 10 mK and IK. We analysed our config- 
urations according to the method described above along 
the symmetry directions [00£] and at a resolu- 

tion of 0.1 reciprocal lattice units. In our analysis, it was 
necessary to account for the classical population statis- 
tics inherent to Monte-Carlo simulations, rather than the 
Bose-Einstein populations implied by Eq. JHJ. The first- 
order finite-size effects given by the zone-centre spin- wave 
frequencies were subtracted from the raw data; these 
were very small for the T = 1 K data, and essentially 
non-existent for the lower temperature set. The error 
associated with each point on the dispersion curves was 
estimated by comparison of the energy values obtained at 
symmetry-equivalent wave- vectors. Our results are com- 
pared with the theoretical dispersion curves in Fig.^ 

The excellent agreement between observed (Monte- 
Carlo generated spin configuration, analysed via the spin 
coordinate matrix) and calculated [Eq. (fTTjf) ] spin-wave 
dispersion is strong evidence for the applicability of the 
theoretical approach described in this paper. In this test 
case, we simply recover the very information that we used 
to drive the Monte-Carlo simulations — namely the value 
of J and the nature of the spin interactions. But the 
key result is that these quantities were extracted solely 
from a set of individual spin orientations. The analy- 
sis is blind to the method with which these were gen- 
erated, and indeed the particular interaction parameters 
employed. We have shown that these parameters may be 
recovered quantitatively from atomistic configurations. 
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IV. DISCUSSION AND CONCLUSIONS 

Our general approach could be automated very easily 
in the form of a computer program, and our limited expe- 
rience of implementing the analysis in this way has shown 
that the calculation of spin- wave dispersion curves of rel- 
atively simple systems can be carried out using a conven- 
tional desktop computer over very manageable (< lh) 
timescales. The generation of sufficiently many configu- 
rations to yield results with acceptable error margins will 
in general demand significantly more extensive computa- 
tional resources. We have found that for simple systems, 
one requires in the vicinity of 500 independent configu- 
rations to give acceptable results. 

When analysing SD simulations, one has access to 
more than spin orientations alone, and the use of this 
additional information can help improve the quality of 
the spin-wave dispersion curves obtained. Given that the 
simulations calculate changes in orientation as a function 
of real time, it is possible to calculate the individual spin 
momenta across a given configuration. Our method of 
analysis may be extended by calculating — in addition to 
the spin displacement matrix S(k) — the related spin mo- 
mentum matrix S(k), constructed from the momentum 
variables (j ± (j£, t) . Then it is easily shown that the ra- 
tio of the eigenvalues e(k, v) of S(k) to the eigenvalues 
e(k, v) of S(k) is related to the spin- wave mode frequen- 
cies: 

9 , e(k, v) 

While computationally more intensive, the calculation 
of spin-wave frequencies in this manner yields results of 



a higher quality, particularly in the form of the long- 
wavelength acoustic modest The origin of this improve- 
ment is particular to the nature of SD simulations, where 
the timescale sampled during the simulation might not 
be sufficient to ensure equipartition of energy between all 
normal modes. The formalism of Eq. (jl 1|> ensures that 
one evaluates the extent of spin displacement in the con- 
text of the population of a particular spin-wave mode, 
rather than the expected population as given by the over- 
all extent of spin-wave excitation. This extension is of 
course irrelevant in the analysis of ergodic simulations, 
where there is no real concept of time. 

In conclusion, we have shown how a relatively straight- 
forward extension of standard spin-wave theory can be 
used to obtain a quantitative link between the energies 
of the various spin-wave modes in a material and the 
individual spin orientations one might observe in "atom- 
istic" models of its behaviour. The method we describe 
has the particular advantage that it allows the extraction 
of spin-dynamical quantities without any prior assump- 
tion of the form of the various spin interactions. This 
model-independence would be of significant advantage 
when using the technique to study materials for which 
little is known a priori about the nature of the magnetic 
interactions, or where one wishes to avoid the presump- 
tion of a particular interaction model. 
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